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Abstract. The molecule solution of an equation related to the lattice Boussinesq equation is 
derived with the help of determinantal identities. It is shown that this equation can for certain 
sequences be used as a numerical convergence acceleration algorithm. Numerical examples with 
applications of this algorithm are presented. 
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1. Introduction. In recent years, it has been found that integrable systems are 
closely connected to certain numerical algorithms. This observation allows a fresh 
look at the research in both fields, and a lot of interesting work has been done lately. 

For example, one step of the QR algorithm is equivalent to the time evolution 
of the finite nonperiodic Toda lattice |3D]. Moreover, Wynn's celebrated £-algorithm 
[55] is nothing but the fully discrete potential KdV equation [HI [53] . The continuous- 
time Toda equation leads to a new algorithm for computing the Laplace transform of 
a given analytic function The discrete Lotka-Volterra system has applications 

in numerical algorithms for computing singular values [161 115[ 132] . As far as the links 
between integrable systems and convergence acceleration algorithms are concerned, 
more results have recently been achieved (see, e.g., [H [TU [T^l [H] ) . 

On the other hand, rapid progress has been made recently in the study of discrete 
integrable systems. As a result, many new examples of discrete integrable systems 
have been found, among them the lattice Gel'fand-Dikii hierarchy [33] [23]. However, 
to the best of our knowledge, nothing has been done so far on designing convergence 
acceleration algorithms via new discrete integrable systems. This is what we want to 
do in this article. Our starting point will be the lattice Boussinesq equation [221 123] , 
and we want to use it for the design of a new convergence acceleration algorithm. 

In view of the importance of both discrete integrable systems and sequence trans- 
formations in the context of convergence acceleration algorithms, it might be helpful 
- before discussing more specific details - to first mention some basic facts about 
discrete integrable systems and about sequence transformations, and to explain how 
discrete integrable systems can be a source of inspiration for sequence transformations. 

Discrete integrable systems could be considered as a specific class of discrete sys- 
tems which possesses an important property of what is usually called integrability. In 
particular, the existence of a Lax pair and a r- function are two of the most impor- 
tant features of integrability shared by some famous numerical algorithms, such as 
Rutishauser's gd-algorithm ^5] and Wynn's e-algorithm [3^. Let us also remark that 
corresponding to different boundary conditions there are different versions of discrete 
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integrable systems available. An example is the famous Toda equation 



(1.1) 



In the infinite chain case with fc = 0, ±1, ±2, . . ., we call (|l.ip the infinite Toda lattice 
equation. Under periodic boundary condition x^+Ko = with fc = 0, ±1, ±2, . . . and 
fixed Kq > 0, (jl.ip is referred to as the periodic Toda lattice equation. If A: = 0, 1, 2, . . . 
with the boundary condition XQ{t) — — oo, we call (jl.ip the semi-infinite Toda equation 
or infinite Toda molecule equation. If XQ{t) — —oo and xjv+i ~ +oo, we call (jl.ll) the 
finite nonperiodic Toda equation or finite molecule Toda equation. 

In connection with convergence acceleration algorithms, we are exclusively in- 
terested in the semi-infinite or infinite molecule case of (jl.ip . corresponding to the 
semi-infinite or infinite Toda molecule equation, and the solutions obtained in this 
way are called molecule solutions. As it turns out, our molecule solutions are closely 
related to sequence transformations. 

We now return to the lattice Boussinesq equation, which is the second equation of 
the lattice Gel'fand-Dikii hierarchy [22l [23| . The multisoliton solutions of the lattice 
Boussinesq equation defined on the elementary square were derived in |13[ I31j . An 
ultradiscrete lattice Boussinesq equation and an alternate form of the discrete poten- 
tial Boussinesq equation had been proposed and the multisoliton solutions of both 
equations had also been obtained [17]. The lattice Boussinesq equation is expressed 
as [221 [23] 



p-q+ - u]i, p-q+ «r+' - "i+r 

+ uTip-q + uT+'-uT+^) 

= (2p + q)iul1, + u^^X') -iP + muT^' + uT+\'), 



which is equivalent to 



p3 _ q3 p3 _ q3 



If we now set 



p-q + uT+t'~ur+, p-q + ur+'-u^lV 

= {p + 2q + <i - - q + _ ^jn+l) 

-ip + 2q + ur- u^lY){p - q + - ^II^i)- (1-2) 

= uj" - pi - qm, p^ -q^ = 1, 



we obtain from (|1.2I) 

1 1 



m+l -m -m+2 _ -m+1 



With the help of the variable transformations 

n= -I, k = m + l, UJ! = u™. 
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we obtain from (|1.3|) 



k+l 



Tjn 

1 



fe+3 



^3) 



n+2 
k 



k+l 



7-rn+l\ 
'-'k+l ) 



1 



^k+2 



Tjn 



u. 



n+2 



'k+2 '-'k+2 

This relationship can be simphfied further yielding the following equation: 

1 



Tjn 



^n+l 



{ur+ 



(1.4) 



'k+2 '^k+2iy^k+l ^k+ll 

In this article, we will first derive the molecule solution of the two-dimensional 
difference equation ()1.4p . Then we will show that the resulting equation can be used 
for the acceleration of the convergence of computationally relevant sequences. 

Our article is organized as follows: In section [21 we will derive the molecule 
solution of (|1.4p with the help of determinantal identities. In section|31 we will provide 
a highly condensed review of the most basic features of sequence transformations. 
In section |4l we will use the results from section [2] to construct a new sequence 
transformation. We will also show that this transformation can be implemented by 
the lattice equation (|1.4I) with given initial values. In section [SJ applications of this 
algorithm are presented. Section [5] is devoted to conclusions and discussions. 

2. Molecule solution of the lattice equation (|1.4p . In this section, we study 
the molecule solution of the lattice equation ()1.4p by Hirota's bilinear method |14j . 
With the help of the dependent variable transformation 



m = 



we obtain the bilinear form of (jl.4p 



F1}G 



•n+l 



pn ■ 
^k 



Fl 



/^n+l T^7i+1 /^n 
h2'~'k-l ^ ^k-1 '-'k+2 



^k ^k+2 



pn+lp 



k+2 



^k+1 ^k > 
Tpn T^n+l 
^k+3^k-l ■ 



(2.1) 
(2.2) 
(2.3) 



Equation (|2.3p can be derived from (j2.ip - p.2p by eliminating the G"s. 
Set 



«'fc(Wn) 



v„ 



Vn+1 



Vn+k-1 
A'^Vn+k-l 
A**W„+fc_i 



A2fe-2„, A2fe-2„, 

A W„ A Vn+1 



^2k-2 



^'_i(w„) = 0, ^'o(w„) = 1, 



and 



^kiVn) = 



$-l(w„) 



n 

Vn 

A2w„ 

A4t;„ 

. 2fe-4, 



n + 1 

Vn+1 
A^Wn+l 
A'^Wn+i 

. 2fe-4, 



A^'^-S„ A-'--Vn+i 

0, $o(Wn) = 1- 



Vn+k-1 



n + k-1 

Vn+k-l 
A'^Vn+k-1 
A'^Vn+k-1 



Vn+k-1 



1,2, 



fc = l,2,.. 
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The solution of an initial value problem related to equations (|2.ip - (|2.3p is given 
below. 



Theorem 1 Given the initial values 



i^i" - F2" = F3" = 1, F^ = AS,„ G^^^O, G^ = n, = Sn- (2.4) 
the solution of the bilinear equations (I2.1I) - (|2.3I) can be expressed as follows: 



= vl/,_i(A35„), 
F3",+i = vI/,(A5„), 
= ^^(A'^r.), 



G^, = vi/,(5„), (2.5) 
G^Vi = -*,_i(A45„), (2.6) 
G^W2 = *fc+i(A^"). (2.7) 



Proof: We consider the case k = 3m in (|2.1I) - (|2.3I) . First, we prove the validity of the 
relationship 



Define 



Di = 



1 

A2 5„ 



3m ^3m 



1 



n+1 



A i^n+l 



1 



n+m 



+m 



(2.8) 



Then we have the relations 



D, = %,,{AS„)=FZ,+^, (2.9) 

i?i(l,2|l,TO+l) -*„_i(A2S„+i) =F^+li, (2.10) 

Dl(l|l) = ^,n{Sn+l) = Gl+\ (2.11) 

i?i(2|m + 1) - *„_i(A35„) = F3"„, (2.12) 

D^{l\m + l) = ^„.{Sn) = G^^, (2.13) 

i?i(2|l) = vI/„,_i(A35„+i) = ^^3"+\ (2.14) 



where Di[j\k) and Di{j, k\p, q) are mth-order and (m — l)th-order determinants ob- 
tained by eliminating the j'th row and the fcth column from the Di and by eliminating 
the jth and fcth rows and the pth and qth columns from the determinant Di , respec- 
tively. 

From the above results, we see that the bilinear equation p.Sp is equivalent to the 
Jacobi identity [10] 

i:)ii:)i(l,2|l,TO+ 1) = Di(l|l)i:>i(2|m + l) - L'i(1|to + 1)Di(2|1). 

Next, we prove the validity of another relationship 

-'^3m+2'-^3?ri-l -'^3m-l'-^3m+2 — ^?,m+l^Zm- \^-^^) 



According to the assumptions of theorem 1, we have 

lu pn ^n+l 



1 



1 

A'S'n+l 



A 



2 m- 1 



Sn A ™ S'n+1 



n + 1 n + 2 

A'^S'n+l A'^S'n+2 



1 

A I? 



A2m— 1 c 

n + m 

^ *-'n+m 



n+2 



^2m-3c 



n-\-m 



1 

A5„+i 

A'^ S'n+1 



1 

A5„+2 

A^S'n+2 



A 2m— 3 c A 2m— 3 c 

^ Jn+1 ^ 'Jn+2 

n n + 1 

ASn A5„+i 

A'^S'n A^Sn+i 



By Schwein's determinantal identity [1] we obtain 



Ihs 



1 

n 
AS„ 
A35„ 



1 

n + 1 

A5'„+i 

A'^ S'n+1 



A^—SSn A2™-3Sn+l 

ASn+1 AS'n+2 
A'^Sn+i A'^S'n+2 



n+2 



We obviously have 



1 

n 

A5n 

A35„ 



1 

n + 1 

ASn+l 
A'^Sn+i 



1 

ASn+r: 

A3S„+, 



^2m-3c 



n -\- m 

m 



A^^-iSn A^'^-iSn+i ••• A^'^-iS, 



1 

n + TO 
A5„+ 



^2m-3c 



n+m 

A^'n-t-m 
n+m 



1 

n + TO 

m 

A3 c 



A2™-3S„ A2™-35„+i ••• A2™-3S, 



n+m 



6 



1 



1 



1 



A 2m — 2 c 



A35„ 
A55„ 



A^S'n+l 
A^S'n+l 



A35, 
A55, 



n+TTi — 2 
n+m — 2 



A2™-i5„ A2™-i5„+i 

= Vl/,„_i(A^5„) = i^3"„,. 

Thus, ;/is = F^'mFs^tli^ and ((2l^ is proved. 
Next, we prove the third relationship: 



A2™-iS'„ 



+m-2 



Define 



1 

A25„ 



A2 5, 
A^S 



n+l 
n+1 



A25, 
A^S 



Then we have the relations 



D2 = ^m{A^Sn) 

L»2(l,m + l|l,m 



n+l 



-^3m+3' 



A2™5„ 



1) = (A^5„+i) = K 



+m 



— 



z?2(i|i) = *™(A25„+i) = ^^3'';„+2, 



3r; 

7^2(™ + 1|™ + 1) = ^-^-iCA^SJ = i^3': 
i^2(l|m + l) = *™(A25„) 
i^2(™ + 1|1) = *™-i(A35„+i) = 



3m+2i 

Tl+1 

3m ■ 



(2.16) 



(2.17) 
(2.18) 
(2.19) 
(2.20) 
(2.21) 
(2.22) 



From the above results, we see that the bilinear equation (|2.16p is nothing but the 
Jacobi identity [TU] 

i:)2i:'2(l,m + l|l,m + 1) = i:»2(l|l)-D2(m + l|m + 1) - D2(l|m + l)D2(m + 1|1). 
The proof of other cases of p.ll) - (|2.3l) can be obtained in a similar way. □ 



3. General properties of sequence transformations. Many calculations 
produce results that are actually sequences whose rate of convergence is governed 
by one or several parameters. Unfortunately, it often happens that the rate of con- 
vergence of such a sequence {Sn} is so slow that the determination of a sufficiently 
accurate approximation to its limit S = Soa by increasing the index n does not work 
in practice. Another frequently occurring problem is that such a sequence {Sn} does 
not necessarily produce a convergent result as n — > oo even if it actually corresponds 
to a meaningful mathematical object with a well-defined numerical value. 



7 



In such a situation, it can be extremely helpful to apply a so-called sequence 
transformation, which transforms the original sequence {Sn\ to a new sequence {5*^} 
with hopefully better numerical properties according to 

T: {5„} ^ {S'J . 

In rudimentary form, sequence transformations have been known for centuries. Their 
older history is reviewed in an article [7 and a monograph [SJ by Brezinski. More 
recent developments are discussed in two articles by Brezinski [5j 6J. There is also an 
extensive bibliography by Brezinski [4, containing more than 6000 references up to 
1991. 

The active research on sequence transformations is documented by the fact that in 
recent years quite a few specialized monographs or longer reviews have been published, 
for example the ones by Brezinski and Redivo Zaglia [9], Sidi [28], and Weniger [33] . 
Numerous other references can be found in |35l Appendix B]. 

Sequence transformations try to achieve an acceleration of convergence or a sum- 
mation in the case of divergence by purely numerical means. Since, however, a com- 
putational algorithm can involve only a finite number of arithmetic operations, a 
sequence transformation T can associate only a finite subset of the input sequence 
{Sn} with an element S'^ of the transformed sequence. 

All the commonly occurring sequence transformations T can be represented by 

(n) 

an infinite set of doubly indexed quantities Tj, with fc,n e Nq. In the literature, 
the superscript n typically indicates the minimal index occurring in the finite sub- 
set {Sn, Sn+i, . • • , Sn+t\ with I = £{k) of sequence elements which are used for the 

(n) 

computation of the transform T^; , and the subscript k - usually called the order of 

(n) (n) 

the transformation - is a measure for the complexity of T^; . Moreover, the are 
gauged in such a way that Tg"'' always corresponds to an untransformed sequence 
element according to Tg"'' = Sn- 

The basic assumption of all sequence transformations is that the elements of a 
slowly convergent or divergent sequence {Sn}, which could be the partial sums 5„ = 
Sfc=o of an infinite series, can for all indices n be partitioned into a (generalized) 
limit S and a remainder or truncation error i?„ according to 5'„ = 5 -I- Rn- If the 
sequence {Sn} converges to its limit S, the remainders i?„ can be made negligible by 
increasing n as much as necessary. But many sequences converge so slowly that this 
does not work in practice. Increasing the index n also does not help in the case of a 
divergent sequence. 

Alternatively, one can try to improve convergence or accomplish a summation by 
computing approximations to the remainders i?„ which are then eliminated from the 
sequence elements Sn, yielding a new sequence with elements S'„ = S + R'^. At least 
conceptually, this is what a sequence transformation tries to do. 

Some transformations - for example the so-called Levin-type transformations dis- 
cussed in [34] - possess closed form expressions. However, the vast majority of all 
known sequence transformations are defined by a recursive scheme. 

The probably best known example of such a transformation is Wynn's e-algorithm 
|36| . which corresponds to the following nonlinear two-dimensional recursive scheme: 
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^(") _ 



-(") _ 

-fc+1 



= e 



0, £, 

(«+l) 



(n) 




1 



fc-1 



.(n+1) 
'-k 



-(») ' 



n e No, 



(3.1a) 
(3.1b) 



The elements £2fe'' with even subscripts provide approximations to the (generahzed) 
hmit S of the sequence {Sn} to be transformed, whereas the elements e2fe+i '^ith odd 
subscripts are only auxiliary quantities which diverge if the whole process converges. 
A compact FORTRAN 77 program for the e-algorithm as well as the underlying com- 
putational algorithm is described in [33l section 4.3], and in [25l p. 213], a translation 
of this FORTRAN 77 program to C can be found. 

The two-dimensional recursive scheme p. II) for the e-algorithm was derived by 
making some assumptions about the mathematical nature of the truncation errors 
or remainders i?„ = Sn — S oi the elements of the sequence {Sn} which is to be 
transformed. Other sequence transformations defined by recursive schemes were also 
constructed by trying to eliminate the remainders of appropriate model sequences. 

But a converse approach is also possible. We can take an equation in two discrete 
variables - for example, one based on integrable systems - and analyze if and under 
which conditions this equation can be used as a starting point for the construction of 
a sequence transformation. 

Let us conclude this section with some useful terminology that is typical of the lit- 
erature on sequence transformations. Assume that a sequence {S",!}, which converges 
to some limit S, satisfies 

Sn+l—S Rn+l /„ „^ 

lim -F^ 7^ = bm — — = p. (3.2) 

n— i-oo i„ — i n— s-oo tin 

If < IpI < 1 holds, we say that the sequence {Sn} converges linearly, if p = 1 holds, 
we say that this sequence converges logarithmically; and if p = holds, we say that it 
converges hyperlinearly. Of course, jpj > 1 implies that the sequence {Sn} diverges. 
Simple examples of linearly, logarithmically, and hyperlinearly convergent sequences 
are the partial sums of the geometric series 1/(1 — z) = Yl^^a'^'^i '^^ Dirichlet 
series C,{s) = X]l^o(^ "I" 1)"" ^'^'^ ^^"^ Riemann zeta function, and of the power series 
exp(z) — X^l^o ^'^ 1^^- foi' the exponential function, respectively. 

Let us assume that two sequences {Sn} and {S'n} converge to the same limit S. 
We say that the sequence {S'n} converges more rapidly than {Sn} if 

lim = 0. (3.3) 

n-i-oo bn — J 



4. A convergence acceleration algorithm. In this section, we propose a new 
sequence transformation and show that this transformation can be derived via the 
lattice equation \\A\ . 

We now consider a new sequence transformation defined as the following ratio of 
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determinants: 



T, 



(n) 



Sn 


Sn+1 

A'^5'^^_i_i 
A''S'„+i 


Sn+k 
A'^ Sn+k 
A'^Sn+k 


A2^5„ 


A'^''Sn+l 


A'^'^Sn+k 


1 


1 

A^S'n+l 
A'^Sn+1 


1 

A^Sn+k 
A'^Sn+k 


A^'^Sn 


A^'^5,1+1 


■ A^'^Sn+k 



(4.1) 



Obviously this transformation is a particular case of the .E-transformation which 
was first derived by Schneider [57] and later rederived independently by Havie [TT] 
and Brezinski From the kernel of the £'-transformation [5], the kernel of the 
transformation (j4.ip is obtained in the theorem below. 

Theorem 2 A necessary and sufficient condition that t/"-* S for all n is that 



Sn = S + fliA Sn + a2A Sn 



-akA^'^Sn, 



where at with i = 1, . . . , fc are some constants. 

Theorem 1 implies that the transformation (j4.ip can be implemented via the 
lattice equation 



Tjn 



1 



with initial conditions 



k 



1,2,.. 



n = 0, 1,.. 



n = 0, 1, 



(4.2) 
(4.3) 



We have 



rpin) _ jjn 
^k — ^3fe+3- 



5. Numerical experiments. In this section, we will show how the convergence 
of some of some slowly convergent example sequences can be accelerated by the new 
algorithm 

Example 1. We consider the linearly convergent sequence 



Sn = 2" sin 



(2") ' 



which converges to 5 = tt = 3.14159 26535 89793 . . .. The corresponding transforma- 
tion results are presented in Table 5.1. 
Example 2. We consider the sequence 



(-1) 



fe-1 



k=l 
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Table 5.1 
Numerical results of example 1. 



n 


^0 






rp(,n) 
^3 


rp{n) 


1 


2.00000 


3.16790 51916 


3.14158 12622 


3.14159 26537 


3.14159 26536 


2 


2.82843 


3.14304 69467 


3.14159 24821 


3.14159 26536 




3 


3.06146 


3.14168 10168 


3.14159 26509 






4 


3.12145 


3.14159 81382 


3.14159 26535 






5 


3.13655 


3.14159 29958 


3.14159 26536 







Table 5.2 
Numerical results of example 2. 



n 


rp{n) 
^0 




rp(n) 


rp(n) 
^3 


rp(n) 




1 


1.00000 


0.70588 


0.69381 


0.69318 36537 


0.69314 92236 


0.69314 72961 


2 


0.50000 


0.68817 


0.69294 


0.69313 74926 


0.69314 66861 


0.69314 71544 


3 


0.83333 


0.69557 


0.69323 


0.69315 03633 


0.69314 73258 


0.69314 71876 


4 


0.58333 


0.69178 


0.69311 


0.69314 59611 


0.69314 71311 




5 


0.78333 


0.69399 


0.69316 


0.69314 77054 


0.69314 71995 





of the partial sums of the ahernating series, which converges to 5* = In 2 = 0.69314 71805 . . .. 
The corresponding transformation results are presented in Table 5.2. 
Example 3. We consider the logarithmically convergent sequence 

n ^ 

fe=i 

of the partial sums of the Dirichlet series for C(2), which converges to S' = ^ = 
1.64493 406684 . . .. The corresponding transformation results are presented in Table 
5.3. 
Note 

The results given above show that our new algorithm accelerates convergence in the 
first two cases but not in the third case. As is well known |9], Wynn's e-algorithm, 
which corresponds to the discrete potential KdV equation, can accelerate linear con- 
vergence and the convergence of alternating series but fails to accelerate logarithmic 
convergence. Apparently, the discrete Boussinesq equation and the discrete potential 
KdV equation have similar properties from the integrable systems' point of view. We 
therefore guess that the algorithm based on (|1.4I) related to the lattice Boussinesq 
equation has the similar acceleration properties as the e-algorithm. The numerical 
results given above are in agreement with our conjecture. 

6. Conclusions and discussions. In this article, we construct the molecule so- 
lution of (jl.4|) related to the lattice Boussinesq equation by Hirota's bilinear method. 
We show that this equation can be used as a numerical convergence acceleration 
algorithm. Our numerical experiments show that this algorithm is effective for lin- 
early convergent sequences and alternating series but fails in the case of logarithmic 
convergence. Also we have studied the confluent form of (|1.4I) in [29]. 

It is known that the fully discrete potential KdV equation and the lattice Boussi- 
nesq equation are the first two equations of the Gel'fand-Dikii hierarchy. Now we 
know that both of the equations have connections with convergence acceleration al- 
gorithms. Therefore, it is a natural idea to investigate further whether higher order 
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Table 5.3 
Numerical results of example 3. 



n 


rp{n) 




rp(n) 


rp(n) 




rp(n) 


rp{n) 


rp{n) 


1 


1.00000 


1.38462 


1.49536 


1.54487 


1.57198 


1.58872 


1.59990 


1.60782 


2 


1.25000 


1.45686 


1.52776 


1.56266 


1.58298 


1.59608 


1.60511 


1.61164 


3 


1.36111 


1.49794 


1.54871 


1.57512 


1.59112 


1.60175 


1.60925 


1.61486 


4 


1.42361 


1.52436 


1.56334 


1.58432 


1.59738 


1.60625 


1.61261 


1.61721 


5 


1.46361 


1.54276 


1.57412 


1.59138 


1.60234 


1.60990 


1.61542 


1.62011 



members of the Gel'fand-Dikii hierarchy have relationships to other convergence ac- 
celeration algorithms. We will consider these problems in the future. 
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